Criticality and isostaticity in fiber networks 
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The rigidity of elastic networks depends sensi- 
tively on their internal connectivity and the na- 
ture of the interactions between constituents. 
Particles interacting via central forces undergo 
a zero-temperature rigidity-percolation transition 
near the isostatic threshold, where the constraints 
and internal degrees of freedom are equal in num- 
beiLUSJ. Fibrous networks, such as those that form 
the cellular cytoskeletorpEl become rigid at a 
lower threshold due to additional bending con- 
straints. However, the degree to which bending 
governs network mechanics remains a subject of 
considerable debateSHSl. We study disordered fi- 
brous networks with variable coordination num- 
ber, both above and below the central-force iso- 
static point. This point controls a broad crossover 
from stretching- to bending-dominated elasticity. 
Strikingly, this crossover exhibits an anomalous 
power-law dependence of the shear modulus on 
both stretching and bending rigidities. At the 
central-force isostatic point — well above the rigid- 
ity threshold — we find divergent strain fluctua- 
tions together with a divergent correlation length 
£, implying a breakdown of continuum elasticity 
in this simple mechanical system on length scales 
less than £. 

There are numerous examples of stiff-fiber networks, 
ranging from carbon nanotube gels at the small scale to 
felt and paper at the macroscopic scale^^. In addition, 
critical biological components such as the mtra-cellular 
cytoskeleton and extra-cellular matrices of collagen and 
fibrin are composed of such network*^. Here, we use 
a combination of simulations and effective medium the- 
ory we study the elasticity of disordered fiber networks 
composed of straight, stiff filaments organized on a trian- 
gular lattice in 2D and face centered cubic (FCC) lattice 
in 3D, as illustrated in Fig. [I] Undiluted, these networks 
have a coordination number z rn — 6 (triangular lattice) 
and z m — 12 (FCC), placing them well above the central- 
force (CF) isostatic point z c = 2d in d dimension d 1 * 17 !. We 
explore the effects of network connectivity — both above 
and below z c — by removing filament segments between 
vertices with a probability 1 — p. 

The mechanical response of the fibers in the network 
is determined by their bending rigidity re and stretch- 
ing modulus /x. For small deformations, the stretching 
and bending energy of the network can be expanded to 
quadratic order in the displacements u^ from the unde- 




FIG. 1: Fiber networks arranged on 2D and 3D lat- 
tices A small section of a sheared diluted triangular network 
near isostaticity with relatively stiff filaments (re = 10 1 in 
units of /j,£o) (a) and floppy filaments (re = 10 -5 ) (b). The 
deviation of the local deformation from a uniform deformation 
is indicated by color, where blue corresponds to a uniform or 
affine deformation and red corresponds to a highly non-affine 
deformation, (c) An example of a small section of the diluted 
FCC network at p = 0.7. To probe the mechanical properties 
of this network we shear the Ill-plane (shown on top) along 
the direction of one of the bond angles in this plane. 

formed reference state at each vertex i, 

-^stretch = o T~ ^ 9ij ( u ij ' ^ij) 2 (1) 

w) 

£bcnd = 9n9jk [(Ujfc - Uij) x f^] 2 , (2) 

(ijk) 

where £q is the lattice spacing, Uy = Uj — U; and 
is the unit vector oriented along the ij-ih bond in the 
undeformed reference state. Here, gij — 1 for present 
bonds and gij — for removed bonds. The summation 
extends over neighboring pairs of vertices in the stretch- 
ing term [Eq. Q], and over coaxial neighboring bonds in 
the bending term [Eq. Q]. Thus, in our networks the 
cross-links at each vertex are freely hinging. 

To investigate the mechanical response of a network, 
we calculate its shear modulus G numerically. The di- 
luted networks exhibit a finite shear modulus well be- 
low the CF isostatic point (expected at p c — 2/3 in 2D 
and p c = 1/2 in 3D), as shown in Fig. [2^,,b; G van- 
ishes at a re-independent rigidity percolation point lo- 
cated at pb = 0.445 ± 0.005 (2D triangular lattice) and 
p b = 0.275 ± 0.005 (3D FCC lattice), consistent with a 
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FIG. 2: Mechanics and non-affine strain fluctuations The shear modulus G in units of fJ,/£o as a function of p for a 
range of filament bending rigidities « for the 2D triangular lattice (a) and the 3D FCC lattice (b). The EMT calculations for 
the 2D triangular lattice are shown as solid lines. The non-affinity measure F is shown as a function of p for various values of 
k for the 2D triangular lattice (c) and the 3D FCC lattice (d). The values of k in units of fiEg are 10° (green), 1CT 3 (cyan), 
10" 4 (red) and 10" 6 (blue). 



floppy mode counting argument that includes the bend- 
ing constraints (suppl. info.). Plots of G versus p for 
different k are shown for the triangular and FCC lattices 
in Figs. [2] a,b. For p > p c , G approaches a nearly k- 
independent stretching dominated limit with G ~ /i. In 
contrast, for p < p c , G falls off reaching a bending domi- 
nated limit with G ~ n as p — > pb- The most interesting 
behavior occurs near p c as a function of k. There is a 
stretch dominated regime at large k and bending domi- 
nated one at small k with a broad intermediate crossover 
regime with G depending on both k and [i. 

To gain insight into the mechanical behavior of our 
models, we developed a new effective medium theory 
(EMT) or coherent potential approximation (CPA^pH*! 
for lattices with bending forces, which we discuss in more 
detail in the methods section, whose results for G for dif- 
ferent k are shown in Fig. [2^,. These results overestimate 
the rigidity percolation point p^. Nonetheless, this model 
captures the essential features of the simulations well, 
including the crossover between stretching and bending 
dominated regimes close to p c . Our EMT theory predicts 
that when /c//z <S! Ap, G can be expressed in the vicinity 
of p c in the scaling form 



G = ^\Ap\fg ± (-\Ap\ 



(3) 



where / = /emt = 1 and < 



i>EMT = 2 are, respectively, 



the rigidity and crossover critical exponents. This scaling 
form is analogous to that for the conductivity of a random 
resistor networlPS with bonds occupied with resistors of 
conductance <t> and a < with respective probabilities p 
and (1— p). When y <C 1, G+{y) ~ const, and 0-(y) ~ y 
implying G ~ p\Ap\f for Ap > and G — K|Ap|^ _<;4 for 
Ap < 0. In the opposite limit (k//i) ^> \Ap\^, G must 
become independent of Ap since it is neither zero nor infi- 
nite at Ap = 0. Equation Jl predicts G ~ n 1 r - u F < , 
which reduces of G ~ k 1 / 2 ^ 1 / 2 in the EMT theory. The 
full EMT results for G along with the scaling form valid 
at k/(j. <C |Ap|^ are shown in Fig. [3^l. 

Our simulation data for both 2D and 3D networks are 
well described by the scaling hypothesis in Eq. pH), con- 
sistent with a second-order transition for k = in both 
cases^l. Fig. [3}j shows the results for both the triangular 
case and FCC case s (ins et). As expected from previ- 
ous simulation worlP^EH, we firici a bending-dominated 
regime at small k and a stretching-dominated regime at 
large k. Consistent with the EMT prediction above, we 
find a previously unexpected intermediate regime with 
G - n x ^~ x where x = f/(j> 0.50 ± 0.01 (2D) and 
f/<p ~ 0.40±0.01 (3D). These results are consistent with 
our exponents obtained above (Table 1, suppl. info.). 
While the extent of this intermediate regime is bounded 
from above by the affine modulus, it can extend to arbi- 
trarily small k > as the system is brought closer to CF 



3 




FIG. 3: Scaling analysis of the mechanics and anomalous elasticity Scaling of the shear modulus in the vicinity of the 
isostatic point with the scaling form G|Ap| - ^ = Q± (re|Ap| - ^), with G in units of n/£o, for the mechanical properties of the 
diluted triangular lattice for the EMT calculations (a) and the simulations (b) for a broad range of filament bending rigidities 
(re in units of filf-,: 1CP 1 black, 10~ 2 magenta, 10~ 3 cyan, 10 -4 red, 10 -5 purple and 10 -6 blue). The asymptotic form of the 
scaling function for low re is shown as a dashed grey line in (a). The scaling for the numerical data on the 3D FCC lattice 
is shown as an inset in b. The EMT exponents are /emt = 1, <^emt = 2. In contrast, for the numerical data we obtain 
/ = 1.4 ± 0.1, = 3.0 ± 0.2 (2D) and / = 1.6 ± 0.2 <j> = 3.6 ± 0.3 (3D). The scaling for the numerical data is performed with 
respect to the isostatic point of the finite system for which we find p c (W) ~ 0.651 (2D, 14 7 =200) and p c (W) w 0.473 (3D, 
W—30). (c) The shear modulus as a function of re close to the isostatic point for the triangular lattice (p — 0.643, blue circles) 
and the FCC lattice [p — 0.47, red squares). At low re there is a bending dominated regime Gbond ~ re, at intermediate re 
there is a regime in which stretching and bending modes couple strongly with G ~ (i 1 ~ x k x , where x = 0.50 ± 0.01 (2D) and 
x « 0.40 ± 0.01 (3D). The EMT calculation for k/(i > | Ap| teMT is shown as a solid blue line. 



isostaticity. 

To investigate the nature of the various mechanical 
regimes, we examine the local deformation field in our 
simulations. Several methods have been proposed to 
qua ntify t he deviation from a uniform (affine) strain 
fieloP^SEH Here we utilize a measure for this non-affinity 
given by 

r = ^E[^"H 2 , (4) 

i 

where u| aff ' ) is the affine displacement of vertex i and 
N is the number of vertices. This quantity varies over 
eight orders of magnitude, indicating non-affine fluctua- 
tions that depend strongly on both re and p, as shown in 
Figs.[2j:,d. For stretch-dominated networks (high re), we 
find a monotonic increase in non-affine fluctuations with 
decreasing p, which appear to diverge at pb- Remarkably, 
for smaller values of re, a second peak in T develops at p c . 
Importantly, the development of this peak coincides with 
the appearance of a crossover between the stretching and 
bending regimes (Figs. [2^-d). 

The critical behavior we observe suggests both a di- 
vergence of the non-affine fluctuations according tcP^ 
r = r_|-|Ap| _A and the existence of a divergent length- 
scale £ = £±|Ap| _l/ near the critical point P c for van- 
ishing re. However, the divergence of £ is limited by 
the system size W, which should suppress the diver- 
gence of r. Consistent with this picture, we find that 
the the location of the cusp in the local fluctuations 



r shift towards higher p with increasing W according 
to p c (W) = p c + bW~ 1 l u , with v = 1.4 ± 0.2 and 
p c = 0.659 ± 0.002 (suppl. info.); these values are con- 
sist with previous reports on generic CF networks^. 
In addition, the amplitude of T increases with system 
size (Fig. |4^i) , in quantitative accord with the expected 
finite-size scaling. Specifically, we find a good collapse 
of the simulation data with T = W X / U F T}± (| AplW 1 ^) 
over a range of system sizes, with \jv — 1.6 ± 0.2 and 
v = 1.4 ± 0.2, as shown in Fig. |4Jd. Similarly, the shear 
modulus exhibits finite-size scaling (Suppl. info.) accord- 
ing toG= W- f /"F Gt± (lAplPF 1 /^), as shown in Fig.Hi. 
We obtain a good collapse of the elasticity data using 
f jv = 0.9 ±0.1, along with v determined from the finite- 
size scaling of T (Fig. [4]b and suppl. info), consistent with 
the value of / obtained from the scaling in Fig. [3] Thus, 
we find a scale-dependence of the shear modulus that is 
consistent with critical behavior governed by the CF iso- 
static critical point. Furthermore, the critical behavior 
in these purely mechanical networks is accompanied by 
shear-induced divergent non-affine fluctuations. These 
results imply a breakdown of continuum elasticity below 
the divergent length-scale £. 

The undiluted triangular and FCC lattices we study 
have an average coordination number greater than 2d 
and thus are above the Maxwell central-force isostatic 
threshold. These networks also consist of infinitely long 
filaments. Cutting bonds as we do introduces both finite 
length polymers, as well as lower connectivity, down to 
the CF threshold and below. Cytoskeletal and extracel- 



4 




FIG. 4: Finite size scaling (a) The non-affinity measure V for the 2D triangular lattice at k = for various systems sizes 
W (25 blue, 50 green, 100 red, 150 cyan and 200 purple). Finite size scaling of the non-affinity measure V according to the 
scaling form F = W x/ "J r r,± (ApW 1/u ) (b) and of the shear modulus with the scaling form G = W~ f/v J r G ,±{\Ap\W 1/v ) (c). 
Here Ap = p - p c , where p c = 0.659 ± 0.002. The exponents we obtain are \/u = 1.6 ± 0.2, v = 1.4 ± 0.2 and f/v = 0.9 ± 0.1. 



lular networks can have z as low as 3 (e.g., in branched 
networks) and as high as 6 (in the case of actin-spectrin 
networks), although they typically have a local connectiv- 
ity z ~ 4, where two filaments are connected by a cross- 
link. As a consequence, the CF isostatic point is expected 
to occur for high molecular weight in 2D. We conjecture 
that there is an analogous crossover behavior for such 
networks, including the anomalous scaling behavior for 
the elasticity. In addition, we expect that our results for 
the crossover behavior will apply to bond-bending models 
on similar lattices to ours^MHI f or rigidity percolation 
and network glasses that include bending forces between 
bonds pairs at each network node. Finally, from the per- 
spective of critical phenomena more generally, the kind 
of crossover behavior we find here is in contrast to most 
thermal systems, where a field or coupling constant leads 
to a crossover from one critical system to another, such 
as from the Heisenberg model to the Ising modePn In 
such systems, there is a continuous evolution of the crit- 
ical point that is governed by the crossover exponent (f>. 
Interestingly, we find no such continuous evolution, but 
rather a discontinuous jump in the critical point p c as 
soon as k becomes nonzero. 



I. METHODS 

Simulations The mechanical response of the network is 
determined in our simulations by applying a shear defor- 
mation with a strain 7. This is realized by translating 
the horizontal boundaries to which the filaments are at- 
tached, after which the internal degrees of freedom are 
relaxed by minimizing the energy using a conjugate gra- 
dient algorithm 2 -. To reduce edge effects in our simu- 
lation, periodic boundary conditions are employed at all 
boundaries. The shear modulus of the network is re- 



lated to the elastic energy through G = yj^yd ^1 for a 
small strain 7, where Vq is the area/volume of a unitcel. 
Here W d is the system size, which in our simulations is 
W 2 « 40000 (2D) and W 3 « 30000 (3D), and we used 
strains no larger than 7 = 0.05. 

EMT The EMT maps the diluted random network to an 
undiluted uniform effective medium (EM) with respective 
stretching modulus fx m and bending rigidity n mi which 
are determined self-consistently as follows. In our the- 
ory, k is as a property of the filament connecting neigh- 
boring sites rather than as a site-associated rigidity that 
connects next-nearest neighbor sites. Following standard 
EM procedures, an arbitrary bond is cither replaced with 
probability p by a bond of stretching modulus jU and 
bending rigidity n or removed with probability 1 — p. 
The phonon Green's function after this replacement is 
calculated as a perturbation with respect to the uniform 
effective medium, treating the replaced bond as a scatter- 
ing potential V on the EM Hamiltonian. The EMT self- 
consistency condition requires that the disorder-averaged 
Green's function equals that of the unperturbed EM, i.e., 
that the average T-matrix arising from the perturbed 
bond vanishes, giving us equations determining fi m and 
K m for given p. 

In the EMT scattering potential V, the stretching term 
is simply proportional to fi — /x m if the bond is occupied 
and — /i m if it is removed. The bending terms must, how- 
ever, be treated differently because replacing a bond gen- 
erates two bending terms, both of which involve second- 
neighbor interactions. This can be understood by con- 
sidering 4 sites ijkl along a filament. If one replaces 
bond jk, two bending terms involving second-neighbors 
ijk and jkl are generated in V. The coefficients of these 
two bending terms can be found by considering a com- 
posite filament connecting ijkl that is composed of rods 
with bending rigidity k s between sites jk and K m between 
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sites ij and jk, respectively, where k s = k if the bond is 
occupied and k s = if it is removed. A direct calcula- 
tion of the minimum bending energy yields the effective 
bending rigidity 

/ 1 1 N- 1 
k c = 2[— + —) , (5) 

and thus the coefficients of the two bending terms involv- 
ing ijk and jkl in V is given by k c — n rn . 
To close the EMT self-consistency equation, 

pT( M ,K) + (l-p)T(0,0)=0, (6) 

where T is the T-matrix constructed from the perturba- 
tion of the scattering potential V, a third- neighbor cou- 
pling X m 

(ijkl) 

must be introduced to the EM and to V accord- 
ingly. Thus the EM is characterized by 3 parameters 
(/U TO , K m , A m ), determined by the self-consistency equa- 
tion We obtained asymptotic solutions to the this 
equation for small K in the vicinity of the CF isostatic 
point, in which fj, m can be written into a scaling form 
same as that of Eq. Q by identifying that the shear 
modulus G — v3/-tm/4) and the scaling function is 

G ± (y) ~ |(±l + v/l + 44j//9) 

(8) 

with A ~ 2.413. 

For k//z <C I Ap|^ EMT , to leading order, the value for 
reduces to 3/x|Ap| for Ap > 0, and yK|Ap| _1 for Ap < 0. 
For k/h > |Ap| feMT we get /x m ~ VA ^ 1 / 2 k 1 / 2 . These 
three scaling regimes correspond to three different slopes 
0,1,1/2 in the G\Ap\~ f vs n\Ap\~<f> plot, as shown in 
Fig.§t. 

Effective medium theories for bond-diluted lattices 
with central-force springs are straightforward because the 
springs reside on an individual bond. In contrast, EMTs 
for lattices with bending forces are less so because bend- 
ing forces reside on two bonds and dilution removes only 
a single bond at a time. Our solution is to treat a given 
bond as a filament segment with bending modulus k s . 
The effective lattice bending modulus for neighboring 



bonds with respective bending moduli Kb and n m is given 
by Eq. |5]). This treatment allows us to unambiguously 
remove one bond at a time. The resultant effective theory 
necessarily includes bend-stretch coupling. A previous 
EMT theory^ treated the bending problem by removing 
two bonds at a time. The result was a theory that lacks 
bend-stretch coupling and predicted separate thresholds 
for the development of non- vanishing /i m and K m , which 
is inconsistent with both the numerical and analytical 
EMT results presented here. 
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FIG. 5: Phase diagram The phase diagram for diluted 
super-isostatic networks. Above the rigidity percolation point 
Zb there are three distinct mechanical regimes: a stretching 
dominated regime with G ~ fi, a bending dominated regime 
with G ~ k and a regime in which bend and stretch modes 
couple with G ~ fi 1 ~ x K, x . Here x is related to the critical 
exponents x = f/<j>. We find here that x = 0.50 ± 0.01 (2D 
triangular lattice) and x = 0.40 ± 0.01 (3D FCC). The me- 
chanical regimes are controlled by the isostatic point z c , which 
acts as a zero-temperature critical point. 
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